;; 批量处理CALIPSO VFM数据并生成地理分布图
; 支持处理指定目录下的所有.hdf文件，并自定义经纬度范围

; 输入目录
input_dir = 'G:\CALIPSO\code\python\su\calipso-vfm-visualization-main\data\'

; 输出目录
output_dir = 'G:\CALIPSO\PRO_CODE\output_geo_wrold_track\\'

; === 自定义经纬度范围 ===
min_lat = 15.0   ; 最小纬度（度）
max_lat = 50.0   ; 最大纬度（度）
min_lon = 90.0   ; 最小经度（度）
max_lon = 125.0  ; 最大经度（度）
; ===

; 图像尺寸（英寸）
image_width = 10.0
image_height = 8.0

; 固定输出分辨率（DPI）
output_dpi = 600

; 确保输出目录存在
IF ~FILE_TEST(output_dir, /DIRECTORY) THEN FILE_MKDIR, output_dir

; 查找所有.hdf文件
file_pattern = input_dir + '*.hdf'
file_list = FILE_SEARCH(file_pattern)
n_files = N_ELEMENTS(file_list)

; 检查是否找到文件
IF n_files EQ 0 THEN BEGIN
  PRINT, '未找到HDF文件！请检查输入目录：', input_dir
  RETURN
ENDIF

PRINT, '找到', n_files, '个HDF文件，开始处理...'

; 批量处理每个文件
FOR i = 0, n_files-1 DO BEGIN
  file_name = file_list[i]
  PRINT, '处理文件:', file_basename(file_name)
  base_name = file_basename(file_name)
  file_core = strmid(base_name, 0, strlen(base_name) - 4)
  
  ; 打开文件
  fid = HDF_SD_START(file_name, /READ)
  IF fid LT 0 THEN BEGIN
    PRINT, '错误：无法打开文件', file_name
    CONTINUE
  ENDIF

  ; 读取数据
  data_name = "Feature_Classification_Flags"
  index = HDF_SD_NAMETOINDEX(fid, data_name)
  sds_id = HDF_SD_SELECT(fid, index)
  HDF_SD_GETDATA, sds_id, fcf
  
  data_name = "Latitude"
  index = HDF_SD_NAMETOINDEX(fid, data_name)
  sds_id = HDF_SD_SELECT(fid, index)
  HDF_SD_GETDATA, sds_id, lat
  
  data_name = "Longitude"
  index = HDF_SD_NAMETOINDEX(fid, data_name)
  sds_id = HDF_SD_SELECT(fid, index)
  HDF_SD_GETDATA, sds_id, lon
  
  HDF_SD_ENDACCESS, sds_id
  HDF_SD_END, fid

  ; === 经纬度范围筛选 ===
  lat_1d = reform(lat)
  lon_1d = reform(lon)
  
  ; 查找符合经纬度范围的索引
  lat_indices = WHERE(lat_1d GE min_lat AND lat_1d LE max_lat, count_lat)
  lon_indices = WHERE(lon_1d GE min_lon AND lon_1d LE max_lon, count_lon)
  
  ; 查找公共索引（兼容IDL 8.5，无ISIN函数）
  common_indices = lon_indices[0:0]  ; 初始化
  count_common = 0
  FOR j = 0, count_lon-1 DO BEGIN
    idx = lon_indices[j]
    IF (WHERE(lat_indices EQ idx, cnt) GT 0) THEN BEGIN
      common_indices = [common_indices, idx]
      count_common = count_common + 1
    ENDIF
  ENDFOR
  IF count_common GT 0 THEN common_indices = common_indices[1:*] ELSE common_indices = []
  
  ; 检查有效数据
  IF count_common LT 2 THEN BEGIN
    PRINT, '警告：经纬度范围', min_lat, '-', max_lat, 'N,', min_lon, '-', max_lon, 'E 内数据点过少，跳过此文件'
    CONTINUE
  ENDIF
  
  ; 子集数据
  lat_subset = lat_1d[common_indices]
  lon_subset = lon_1d[common_indices]
  data = fcf(1256, common_indices)  ; 2500m高度数据
  ; ===

  ; 特征类型筛选
  data = data AND '0007'X
  
  ; 数据格式转换
  lat = reform(lat_subset)
  lon = reform(lon_subset)
  data = reform(data)
  ; 修改后的代码（添加自动换行）
;  longname = 'Feature Type at altitude=2500m' + $
;    ' (Lat: ' + strtrim(STRING(min_lat),1) + '_' + strtrim(STRING(max_lat),1) + 'N,' + $
;    ' Lon: ' +  strtrim(STRING(min_lon),1) + '_' +  strtrim(STRING(max_lon),1) + 'E)'
longname = 'Feature Type at altitude=2500m'
  ; 颜色表定义
  levels = 8
  ct = COLORTABLE([[150, 0,   255,  0,   255,  200, 100, 50], $
                  [150,  0,   255,  255,  0,   100, 50,  25], $
                  [150, 255,   0,    0,   0,   255, 255, 125]], $
                  NCOLORS = levels, /TRANSPOSE)

; 生成绘图
; 使用 MAP 函数设置地图区域
; m = MAP('Geographic', LIMIT=[min_lat, max_lat, min_lon, max_lon],TITLE='亚洲地区地图', FONT_SIZE=10, /BUFFER)
;MAP——LIMIT
;Set this property to a four-element vector of the form
;[Latmin, Lonmin, Latmax, Lonmax]
;that specifies the boundaries of the region to be mapped. (Lonmin, Latmin) and (Lonmax, Latmax) are 
;the longitudes and latitudes of two points diagonal from each other on the region's boundary.
;;; 左下的经纬度与右上的经纬度去限制MAP函数
;;; 绘图范围
  m_min_lat = 15.0   ; 最小纬度（度）
  m_max_lat = 50.0   ; 最大纬度（度）
  m_min_lon = 95.0   ; 最小经度（度）
  m_max_lon = 130.0  ; 最大经度（度）
  m = MAP('Geographic', LIMIT=[m_min_lat, m_min_lon, m_max_lat, m_max_lon], FONT_SIZE=9, BACKGROUND_COLOR='light blue', /BUFFER) ; cyan   light blue
;  m = MAP('Geographic', TITLE=file_core, FONT_SIZE=9, /BUFFER)
  t1 = TEXT(0.3, 0.05, longname)   ; longname在图中的位置
  
; === 添加海岸线 ===
; 方法1：绘制海岸线轮廓（推荐）
mc = MAPCONTINENTS(COMBINE=0, $        
                     COLOR='yellow', $  
                     LINESTYLE='solid', $       
                     Thick=2.0)                 
                   ; 添加到地图对象
; 可选：填充大陆区域
mc = MAPCONTINENTS(FILL_COLOR="light gray", TRANSPARENCY=50)  ; dark gray   light green
                                                 
  data = BYTE(data)
  c1 = SCATTERPLOT(lon, lat, OVERPLOT=m, $
                   MAGNITUDE=data, $
                   RGB_TABLE=ct, $
                   POSITION=[0.1, 0.1, 0.83, 0.9],$
                   /SYM_FILLED, SYMBOL='o', SYM_SIZE=0.1)
  mc = MAPCONTINENTS()

  ; 颜色条配置
  cb = COLORBAR(RGB_TABLE=ct, RANGE=[0,8], ORIENTATION=1, BORDER=1,$
                TICKVALUES=[0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5], $
                TICKNAME=['invalid', 'clear', 'cloud', 'aerosol', $
                'strato.', 'surface', 'subsurf.', 'no signal'], $
                TEXTPOS=1, POSITION=[0.85,0.2,0.87,0.8], TITLE=units)

  ; 生成输出文件名（包含经纬度范围）
;  output_file = output_dir + 'geo_' + strtrim(STRING(min_lat),1) + '-' + strtrim(STRING(max_lat),1) + 'N_' + strtrim(STRING(min_lon),1) + '-' + strtrim(STRING(max_lon),1) + 'E_' + file_core + '.png'
  output_file = output_dir + 'geo_wrold_track_'+ file_core + '.png'

  ; 计算像素尺寸（600 DPI）
  pixel_width =  ROUND(image_width * output_dpi)
  pixel_height = ROUND(image_height * output_dpi)
  
  ; 保存图像
  c1.save, output_file, HEIGHT=pixel_height, WIDTH=pixel_width
  PRINT, '已保存:', output_file, ' (', pixel_width, 'x', pixel_height, '像素, 600 DPI)'
ENDFOR

PRINT, '批量处理完成！'
end